Установка R на компьютер состоит из двух частей:
Более подробную инструкцию для разных ОС можно почитать тут.
Обозначения:
\(X_i\) – независимые переменные (covariates)
\(T_i\) – бинарная переменная воздействия (treatment variable): \[\begin{equation*} T_i = \begin{cases} 1, &\text{воздействие на объект i оказано}\\ 0, &\text{воздействие на объект i не оказано} \end{cases} \end{equation*}\]
\(Y_{i1}\), \(Y_{i0}\) – потенциальные исходы (potential outcomes)
Наблюдаемые исходы \(Y_i\) отличаются от потенциальных исходов. Потенциальные исходы являются гипотетическими случайными величинами, когда наблюдаемые исходы являются фактическими случайными величинами. Наблюдаемые исходы являются функцией от потенциальных исходов:
\(Y_i = T_i \cdot Y_{i1} + (1-T_i) \cdot Y_{i0}\)
Действительно при \(T=1\) мы получим \(Y_i = Y_{i1}\), а при \(T=0\) получим \(Y_i = Y_{i0}\).
Тогда наш эффект воздействия для конкретного наблюдения равен разнице между двумя состояниями мира для этого наблюдения (потенциальными исходами):
\(\tau_i = Y_{i1} - Y_{i0}\)
Фундаментальная проблема причинного вывода (Fundamental problem of causal inference):
Чтобы оценить эффект воздействия для конкретного индивида, мы должны знать потенциальные исходы сразу для двух его состояний мира, но реально мы наблюдаем только одно из них – либо \(Y_{i1}\), если индивид подвергся воздействию, либо \(Y_{i0}\), если он ему не подвергался. Оценка индивидуального эффекта требует доступа к данным, которых у нас физически не может быть.
Если с распределением индивидуального эффекта воздействия (treatment effect) работать не получается, будем довольствоваться средними величинами. Например, попробуем рассчитать средний эффект воздействия (average treatment effect):
\(ATE = \mathbb{E}[\tau_i] = \mathbb{E}[Y_{i1} - Y_{i0}] = \mathbb{E}[Y_{i1}] - \mathbb{E}[Y_{i0}] \xrightarrow{p} \frac{1}{N_1}\sum \limits_{i=1}^{N_1}Y_{i1} - \frac{1}{N_0}\sum \limits_{i=1}^{N_0}Y_{i0}\)
Следующий шаг, который мы предпримем, попробуем рассчиать величину эффекта только для тритмент группы – средний эффект воздействия на задействованных (average treatment effect for the treatment group):
\(ATT = \mathbb{E}[\tau_i|T_i=1] = \mathbb{E}[Y_{i1} - Y_{i0}|T_i=1] = \mathbb{E}[Y_{i1}|T_i=1] - \mathbb{E}[Y_{i0}|T_i=1]\)
А теперь то же самое, но для контрольной группы – средний эффект воздействия на незадействованных (average treatment on the non-treated):
\(ATnT = \mathbb{E}[\tau_i|T_i=0] = \mathbb{E}[Y_{i1} - Y_{i0}|T_i=0] = \mathbb{E}[Y_{i1}|T_i=0] - \mathbb{E}[Y_{i0}|T_i=0]\)
Обратите внимание, как и в случае с определением эффекта воздействия на индивидуальном уровне, различные модификации средних эффектов воздействия снова требуют от нас знания обоих потенциальных исходов для каждого наблюдения. Таким образом, и средние, и индивидуальный эффект воздействия нельзя напрямую рассчиать, но мы будем пробовать их оценить.
Самая простая идея для оценки ATE, которая всем придет в голову, взять простую разницу в средних: \(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0]\).
Но тут всё не так просто, после небольших преобразований, с которыми можно ознакомиться в части 4.1.3 учебника мы получим следующее:
\(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0] = \underbrace{\mathbb{E}[Y_1] - \mathbb{E}[Y_0]}_{\text{ATE}} + \underbrace{\mathbb{E}[Y_0|T=1] - \mathbb{E}[Y_0|T=0]}_{\text{Selection Bias}} + \underbrace{(1-\pi)(ATT - ATnT)}_{\text{Heterogeneous treatment effect bias}}\)
Далее мы рассмотрим предпосылки, которые позволяют нивелировать влияние этих двух смещений и получить оценку ATE.
Экзогенность воздействия (Independence assumption) означает, что распределение объекта в тритмент или контрольную группы осуществляется случайно и независимо от его изначальных характеристик. Данная предпосылка обычно обозначается следующим образом \((T_1, Y_0, X)_i \perp T_i\)
Технически для нас это значит следующее:
То есть хорошая рандомизация, а следовательно, и выполнение предпосылок, позволяет нам очистить эффект воздействия от двух типов смещения, в этом случае:
\(ATE = \mathbb{E}[Y_1] - \mathbb{E}[Y_0] = \mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0] \xrightarrow{p} \frac{1}{N_1}\sum \limits_{i=1}^{N_1}Y_{i1} - \frac{1}{N_0}\sum \limits_{i=1}^{N_0}Y_{i0}\)
Эта предпосылка подразумевает выполнение двух вещей. Во-первых, воздействие оказывается только на один объект и внешние эффекты у него отсутствуют. Во-вторых, воздействие гомогенно – существует только один тип тритмента.
В рамках курса нам неоднократно придется прибегнуть к генерации случайных чисел из какого-нибудь закона распределения, как правило, из нормального. Для работы со случайными числами из нормального распределения в R используется несколько функций: dnorm() для плотности вероятности, pnorm() для функции распределения, qnorm() для квантилей распределения и rnorm() для генерации случайных чисел. Почитать подробнее про них и посмотреть примеры можно тут или тут.
Если вы не уверены в синтаксисе какой-то функции, вы можете воспользоваться встроенной справкой в RStudio (справочная информация появится во вкладке Help):
?rnorm
Если вам когда-нибудь понадобятся функции для других распределений, вы можете, конечно, их просто загуглить или ввести в окне Help запрос “Distributions”:
?Distributions
Теперь, когда мы точно знаем, как устроены аргументы у функции rnorm(n;mean;sd), давайте попробуем достать 5 чисел из стандартного нормального распределения \(\mathcal{N} \sim \left(0;1\right)\):
rnorm(5,0,1)
## [1] 0.25331851 -0.02854676 -0.04287046 1.36860228 -0.22577099
Если вы повторяли за мной, то у вас, вероятно, получились другие значения. А теперь я тоже попробую ещё раз:
rnorm(5,0,1)
## [1] 1.5164706 -1.5487528 0.5846137 0.1238542 0.2159416
Как и ожидалось, у меня тоже не совпало с предыдущим результатом. Существует по меньшей мере две причины, почему нам хотелось бы получать одинаковые результаты:
В R существует несколько алгоритмов, позволящих генерировать случайные числа (Random Number Generator (RNG), подробнее – ?RNGkind). По умолчанию используется Mersenne twister (Вихрь Мерсенна), поскольку он работает наиболее быстро. Однако все эти алгоритмы на самом деле детерминированы, поэтому генерируемые ими числа корректнее называть “псевдослучайными”.
Генератор псевдослучайных чисел начинает свою работу с определенной точки в пространстве возможных чисел. Эту точку приянто называть начальное число или seed. Начальное число – это число (или вектор), используемое для инициализации генератора псевдослучайных чисел. Если вы хотите получить одинаковые результаты с помощью генератора случайных чисел, важно установить начальное число. Для этого в R используется функция set.seed(). По умолчанию начальное число не установлено, если ничего не указать, то создается новое из текущего времени и идентификатора процесса на вашем устройстве. Следовательно, по умолчанию разные сеансы дают разные результаты моделирования.
Давайте снова попробуем сгенерировать 5 чисел из случайного нормального распределения, но в этот раз зафиксируем seed:
set.seed(123)
rnorm(5,0,1)
## [1] -0.56047565 -0.23017749 1.55870831 0.07050839 0.12928774
И еще разок:
set.seed(123)
rnorm(5,0,1)
## [1] -0.56047565 -0.23017749 1.55870831 0.07050839 0.12928774
Получилось! В качестве seed можно использовать любое число, главное, чтобы при репликации оно каждый раз было одинаковым.
Почитать про set.seed можно тут на русском и тут на английском.
Симуляции – это “игрушечные” примеры, которые позволят нам на протяжении курса разбираться с разными методами. Игрушечные они потому, что за ними не стоят реальные данные. Данные для симуляций мы будем специальным образом заранее моделировать. Это бывает удобно, когда идеально подходящих данных нет, или они не лежат в открытом доступе. К тому же, живые данные часто могут быть зашумлены из-за других факторов, на которые нам не всегда будет удобно отвлекаться на занятиях. Но с живыми данными мы тоже обязательно будем работать! В том числе, они будут доступны вам в домашних заданиях :)
Теперь смоделируем небольшую симуляцию по мотивам изученого материала.
Давайте смоделируем гипотетическую ситуацию. Мы хотим оценить величину эффекта от использование сайта с расписанием cacs.ws на свободное время студента.
Предположим, что наша экспериментальная выборка состоит из 1000 человек:
N <- 1000
Для простоты у них будет всего две характеристики (\(X\) – возраст и \(Z\) – время в пути от дома до ЭФ), имеющих влияние на потенциальный исход (\(Y_0\) и \(Y_1\) – свободное время). При этом предположим, что реальный эффект воздействия равен \(\tau=15\) минутам. Будем считать, что реальная зависимость потенциального исхода от ковариатов и тритмента выглядит следующим образом:
Сгенерируем \(X\) и \(Z\) случайным образом:
set.seed(123)
X <- runif(N, 18, 25)
Z <- rnorm(N, 60, 20)
# Посмотрим на наши данные, нарисуем сразу 2 графика -- диаграмму разброса и гистограмму
par(mfrow=c(2,2)) # читаем буквально как матрицы -- 2 строки и 2 столбца
plot(X)
hist(X)
plot(Z)
hist(Z)
Также создадим случайные ошибки для каждого из типов индивидов:
set.seed(123)
e0 <- rnorm(N/2, 0, 1)
e1 <- rnorm(N/2, 0, 1)
plot(e0, col="red")
points(e1, col='blue')
Вспомогательно пронумеруем наши наблюдения:
number <- 1:N
И соберем всё в общий датасет:
data <- data.frame(number = number, X=X, Z=Z) # сшиваем столбики в data frame
head(data, 10) # смотрим только первые несколько строк датасета
Что касается рандомизации, то наша глобальная задача, используя разные функции придумать способ, как случайным образом разбить выборку на 2 группы. Существует огромное множество вариантов, как это можно сделать. Мы обсудим лишь несколько из них. На практике вы можете пользоваться как этими вариантами, так и придумать что-то своё.
Первый вариант – просто взять половину наблюдений по их порядку расположения в датасете. Пусть воздействие будет оказано только на первую половину нашей выборки.
T1 <- c(rep(1,N/2), rep(0,N/2))
mean(T1) # проверяем долю тритмента, должна получиться ровно половина
## [1] 0.5
data$T1 <- T1 # добавляем переменную в датасет
Тут мы воспользовались функцией c(...), которая позволяет объединить величины внутри неё в вектор, и функцией rep(x;times) которая повторяет величину \(x\) столько раз, сколько указано в величине \(times\).
Второй вариант, который практически дублирует первый, но немного отличается – мы можем взять каждое второе наблюдение:
T2 <- rep(0:1, N/2)
mean(T2) # проверяем долю тритмента, должна получиться ровно половина
## [1] 0.5
data$T2 <- T2 # добавляем переменную в датасет
Логично, что если мы хотим случайным образом присвоить наблюдениям разные группы, то мы можем начать с того, что присвоим им случайные числа, которые мы потом разными способами переведем в бинарный формат, который уже будет соотвествовать конкретной группе. Протестируем этот способ на примере трех распределений – нормального, равномерного и биномиального.
set.seed(123)
V1 <- runif(N) # генерим вспомогательную переменную из равномерного распределения
T3 <- as.numeric(V1 < median(V1)) #генерим тритмент, отсекая половину выборки по медиане; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
mean(T3) # проверяем долю тритмента, должна получиться ровно половина
## [1] 0.5
data$T3 <- T3 # добавляем переменную в датасет
set.seed(123)
V2 <- rnorm(N) # генерим вспомогательную переменную из нормального распределения
T4 <- as.numeric(V2>0) #генерим тритмент, отсекая половину выборки по знаку вспомогательной переменной; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
summary(T4) # проверяем долю тритмента, должна получиться примерно половина
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 0.505 1.000 1.000
data$T4 <- T4 # добавляем переменную в датасет
set.seed(123)
T5 <- rbinom(N, 1, 0.5) #генерим тритмент, отсекая половину выборки по знаку вспомогательной переменной; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
summary(T5) # проверяем долю тритмента, должна получиться примерно половина
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.493 1.000 1.000
data$T5 <- T5 # добавляем переменную в датасет
experimentlibrary('experiment')
set.seed(123)
rand <- randomize(data, group= c("Treat", "Control"))
T6 <- as.numeric(rand$treatment == "Treat") # преобразовываем переменную в численный формат
data$T6 <- T6 # добавляем переменную в датасет
head(data, 10)
radomizrlibrary('randomizr')
set.seed(123)
T7 <- complete_ra(N=N, m=N/2)
data$T7 <- T7 # добавляем переменную в датасет
head(data, 10)
Существует ряд ситуаций, когда подходы к рандомизации, которые мы разобрали выше, работают не очень хорошо. Например, в вашу выборку добавилось несколько наблюдений, которые по какой-то причине добавились не в конец вашего датасета. В этом случае рандомизация не будет воспроизводимой, а тритмент и контрольная группы будут различаться. Однако есть способ решить эту сложность.
Хэш-функция (hash function) – функция, осуществляющая преобразование массива входных данных произвольной длины в выходную битовую строку установленной длины, выполняемое определённым алгоритмом. Преобразование, производимое хэш-функцией, называется хэшированием.
Хэш-функции применяются в следующих случаях:
Для нас важно то, что с помощью хэш-функции мы сможем взаимно однозначно переводить данные формата строки в число.
Пример хэширования с помощью пакета digest:
library("digest")
hash <- digest('econometrics', algo="murmur32")
hash
## [1] "24b6de5a"
Пакет digest применяет хеш-функцию к произвольным объектам R. В пакете реализовано много разных алгоритмов преобразований, однако мы будем использовать алгоритм “murmur32” (подробнее почитать можно тут). Этот алгоритм совершенно не подходит для криптографических целей, но зато идеально подходит нам, поскольку он 32-битный, что позволяет нам перевести полученный хэш-код в целое число.
Прежде чем начать разбирать пример, немного отвлечемся на техническую полезную вещь – функцию sapply().
Про семейство apply функций рекомендую почитать подробнее на русском тут или тут, а на английском тут.
Прелесть функции sapply() состоит в том, что она позволяет нам избежать громоздких циклов при написании кода и ускорить вычисления благодаря “векторной ориентированности” языка R. Например, функция позволяет нам найти минимальное и максимальное значение для каждой из ковариат:
sapply(list(X,Z), min)
## [1] 18.003257 3.804506
sapply(list(X,Z), max)
## [1] 24.99583 127.80742
То есть функция проводит однотипную операцию (по сути вложенную функцию) над каждым элементом списка. Если кто-то хорошо владеет python, то аналогичной функцией там является map().
Вернемся к нашей симуляции. Изначально мы это все затеяли, чтобы сделать хорошую рандомизацию. Сначала хэшируем номера наших наблюдений. Так делать не очень хорошо, обычно, для хэширования данных используют ФИО и/или СНИЛС, но мы не будем ради этого громоздить генерацию еще одной переменной. Для иллюстрации просто воспользуемся номерами, это тоже сработает.
hashes <- sapply(data$number, function(x) {digest(x, algo="murmur32")}) # хэшируем строчки
hashes[1:20]
## [1] "5e6216f3" "8b8b3789" "ad1bf356" "b2df92b2" "1d1b8e3a" "98d09cab"
## [7] "2e366c5a" "4d110361" "91d6bfe0" "f568e543" "08bd74cd" "e979af4b"
## [13] "a82e5827" "485cbd08" "bc4e5bd8" "95f68430" "7180d744" "da414596"
## [19] "70aa2885" "afd0ffd2"
Далее нужно перевести получившиеся строчки в цифровой числовой. Функция strtoi() конвертирует строковое представление числа, которое хранится в строке, в длинное целое.
result <- strtoi(substring(hashes, 2), base=16)
result[1:20]
## [1] 241309427 193673097 219935574 48206514 219909690 147889323 238447706
## [8] 219218785 30851040 90760515 146633933 158969675 137254951 140295432
## [15] 206461912 100041776 25220932 172049814 11151493 265355218
Далее, используя эти числа, нужно как-то разбить выборку на тритмент и контроль. Будем смотреть на последнюю цифру, если она меньше 5, то наблюдение попадет в тритмент, если больше – в контроль.
T8 <- result %% 10 < 5 # берем остаток от деления на 10 (получится последняя цифра) и сравниваем его с 5
data$T8 <- as.numeric(T8) # переводим из логического типа в числовой и добавляем в датасет
data$T8[1:20]
## [1] 0 0 1 1 1 1 0 0 1 0 1 0 1 1 1 0 1 1 1 0
Проверим нашу рандомизацию на сбалансированность:
summary(data)
## number X Z T1 T2
## Min. : 1.0 Min. :18.00 Min. : 3.805 Min. :0.0 Min. :0.0
## 1st Qu.: 250.8 1st Qu.:19.78 1st Qu.: 46.232 1st Qu.:0.0 1st Qu.:0.0
## Median : 500.5 Median :21.43 Median : 60.576 Median :0.5 Median :0.5
## Mean : 500.5 Mean :21.48 Mean : 60.239 Mean :0.5 Mean :0.5
## 3rd Qu.: 750.2 3rd Qu.:23.23 3rd Qu.: 73.101 3rd Qu.:1.0 3rd Qu.:1.0
## Max. :1000.0 Max. :25.00 Max. :127.807 Max. :1.0 Max. :1.0
## T3 T4 T5 T6 T7
## Min. :0.0 Min. :0.000 Min. :0.000 Min. :0.0 Min. :0.0
## 1st Qu.:0.0 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.0 1st Qu.:0.0
## Median :0.5 Median :1.000 Median :0.000 Median :0.5 Median :0.5
## Mean :0.5 Mean :0.505 Mean :0.493 Mean :0.5 Mean :0.5
## 3rd Qu.:1.0 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.0 3rd Qu.:1.0
## Max. :1.0 Max. :1.000 Max. :1.000 Max. :1.0 Max. :1.0
## T8
## Min. :0.000
## 1st Qu.:0.000
## Median :1.000
## Mean :0.508
## 3rd Qu.:1.000
## Max. :1.000
Посчитаем потенциальные исходы:
Y0 <- 240 - 3*X - Z + 15*0 + e0 # T=0
Y1 <- 240 - 3*X - Z + 15*1 + e1 # T=1
data$Y0 <- Y0 # добавляем переменные в датасет
data$Y1 <- Y1
А также наблюдаемые исходы:
data$Y <- T1*Y1 + (1-T1)*Y0
Когда данные подготовлены, мы забываем, что что-то о них знали :)
С данного момента мы действуем согласно предпосылке, что мы знаем только \(Y\), \(X\) и \(T\).
Поскольку мы знаем, что тритмент распределялся среди наблюдений случайно, мы можем сделать вывод, что предпосылка о независимости воздействия выполнена. Также мы знаем, что данные устроены так, что SUTVA тоже выполнена. Следовательно, можем использовать обычную разницу в средних, чтобы оценить эффнект воздействия:
\(\widehat{ATE} = \overline{Y_1} - \overline{Y_0}\)
ATE_hat <- mean(data$Y[T1==1]) - mean(data$Y[T1==0])
ATE_hat
## [1] 15.61754
Из курса ЭКМ-2 помним, что то же самое можно было бы получить, с помощью обычного парного МНК:
model1 <- lm(Y ~ T1, data=data)
summary(model1)
##
## Call:
## lm(formula = Y ~ T1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -76.984 -14.469 -0.314 13.235 60.998
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 115.0259 0.9094 126.49 <2e-16 ***
## T1 15.6175 1.2861 12.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.33 on 998 degrees of freedom
## Multiple R-squared: 0.1287, Adjusted R-squared: 0.1279
## F-statistic: 147.5 on 1 and 998 DF, p-value: < 2.2e-16
Если обе предпосылки выполнены, то оценка должна быть несмещенной даже без контрольных переменных (ковариат). Попробуем их добавить и сравним оценки эффекта:
model2 <- lm(Y ~ T1 + X + Z, data=data)
summary(model2)
##
## Call:
## lm(formula = Y ~ T1 + X + Z, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9500 -0.4804 -0.0218 0.4296 3.2849
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 238.445186 0.296642 803.8 <2e-16 ***
## T1 14.978356 0.052761 283.9 <2e-16 ***
## X -3.001684 0.013116 -228.9 <2e-16 ***
## Z -0.973141 0.001318 -738.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8341 on 996 degrees of freedom
## Multiple R-squared: 0.9985, Adjusted R-squared: 0.9985
## F-statistic: 2.266e+05 on 3 and 996 DF, p-value: < 2.2e-16
Видим, что результаты устойчивы, величина эффекта практически не изменилась, но наша первая оценка все-таки была слегка завышеная, значит небольшое смещение было. Проверим контрольную и тритмент группу на “похожесть”, то есть проведем баланс ковариатов.
Далее в качестве упраженения посчитаем средние эффекты на разных группах:
ATT <- mean(data$Y1[data$T1 == 1] - data$Y0[data$T1 == 1])
ATT
## [1] 14.96307
ATnT <- mean(data$Y1[data$T1 == 0] - data$Y0[data$T1 == 0])
ATnT
## [1] 14.96307
Получили, то что и должны были получить – ATT = ATnT. Если мы вернемся в раздел, где мы подготовливали данные, то увидим, что в двух потенциальных исходах мы “зашифровали” одинаковый эффект:
Содержательно это иллюстрирует то, что эффект воздействия гомогенен, то есть \(\text{Heterogeneous treatment effect bias} = (1-\pi)(ATT - ATnT) = 0\).
Теперь проверим насколько наша рандомизация оказалась хорошей. Сравним средние значения ковариат в двух группах.
mean(data$X[T1==1])-mean(data$X[T1==0])
## [1] -0.02791759
mean(data$Z[T1==1])-mean(data$Z[T1==0])
## [1] -0.5707093
Это не очень информативный способ, но если мы вспомним, что порядок значений ковариат измерялся десятками, то можно сделать вывод, что группы в среднем похожи.
А теперь проведем тест на разницу в средних с помощью формального t-теста Стьюдента:
t.test(data$X[T1==1], data$X[T1==0])
##
## Welch Two Sample t-test
##
## data: data$X[T1 == 1] and data$X[T1 == 0]
## t = -0.21924, df = 997.51, p-value = 0.8265
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.2777933 0.2219581
## sample estimates:
## mean of x mean of y
## 21.46699 21.49490
t.test(data$Z[T1==1], data$Z[T1==0])
##
## Welch Two Sample t-test
##
## data: data$Z[T1 == 1] and data$Z[T1 == 0]
## t = -0.45032, df = 997.68, p-value = 0.6526
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.057656 1.916238
## sample estimates:
## mean of x mean of y
## 59.95331 60.52402
Согласно значению p-value мы принимаем \(H_0\) и делаем вывод, что группы одинаковые.
library("tableone")
table1 <- CreateTableOne(vars=c("X", "Z"), strata="T1", data=data, test=TRUE)
print(table1)
## Stratified by T1
## 0 1 p test
## n 500 500
## X (mean (SD)) 21.49 (2.04) 21.47 (1.99) 0.827
## Z (mean (SD)) 60.52 (19.86) 59.95 (20.22) 0.653
Можно ли было бы сделать все то же самое за меньшее число строчек кода и время? Конечно, можно. Более того, если бы похожую процедуру вам пришлось бы по каким-то причинам проводить несколько раз, например, когда у вас несколько экспериментов, функция для вас просто была бы незаменимым помощником. Попробуем переписать все то же самое в несколько функций.
generate_data <- function(N){
N=N
X <- runif(N, 18, 25)
Z <- rnorm(N, 60, 20)
e0 <- rnorm(N/2, 0, 1)
e1 <- rnorm(N/2, 0, 1)
number <- 1:N
Y0 <- 240 - 3*X - Z + 15*0 + e0
Y1 <- 240 - 3*X - Z + 15*1 + e1
data <- data.frame(number = number, X=X, Z=Z, Y0=Y0, Y1=Y1)
}
randomization <- function(data, type){
N <- nrow(data)
if (type=='in order'){
T <- c(rep(1,N/2), rep(0,N/2))
data$T <- T
} else if (type=='by turns'){
T <- rep(0:1, N/2)
data$T <- T
} else if (type=='unif'){
V <- runif(N)
T <- as.numeric(V1 < median(V1))
data$T <- T
} else if (type=='norm'){
V <- rnorm(N)
T <- as.numeric(V2>0)
data$T <- T
} else if (type=='binom'){
T <- rbinom(N, 1, 0.5)
}
return(data)
}
outcome <- function(data){
observed_data <- data
observed_data$Y <- data$T*observed_data$Y1 + (1 - data$T)*observed_data$Y0
return(observed_data)
}
estimate <- function(data, response_variable) {
m <- mean(data[data$T==1, response_variable])-mean(data[data$T==0, response_variable])
return(m)
}
balance <- function(data, covariates, type){
if (type=='mean'){
for (i in 1:length(covariates)){
cov <- covariates[i]
m <- mean(data[data$T==1, cov])-mean(data[data$T==0, cov])
print(cov)
print(m)
}
} else if (type=='t.test'){
for (i in 1:length(covariates)){
cov <- covariates[i]
t <- t.test(data[data$T==1, cov], data[data$T==1, cov])
print(cov)
print(t)
}
}
}
А теперь посмотрим, сколько строчек займет наша симуляция с использованием функций:
set.seed(123)
data1 <- generate_data(1000)
set.seed(123)
data1 <- randomization(data=data1, type='in order')
data1 <- outcome(data=data1)
tail(data1,10)
estimate(data=data1, response_variable = 'Y')
## [1] 15.53169
balance(data=data1, covariates = c('X','Z'), type = 'mean')
## [1] "X"
## [1] -0.02791759
## [1] "Z"
## [1] -0.5707093
Ваш исследовательский вопрос может быть таким, что вам интересно оценить воздействия разных типов тритмента, то есть у вас есть несколько экспериментальных групп и одна контрольная. При такой постановке мы хотим проверить не одну, а сразу много статистических гипотез о различии в группах. При проверке любой гипотезы существует вероятность совершить ошибку первого рода (отклонить нулевую гипотезу, если она верна = обнаружить эффект, которого нет). Особенность множественного тестирования гипотез состоит в том, что чем больше гипотез мы проверяем на одних и тех же данных, тем больше будет вероятность допустить как минимум одну ошибку первого рода – эффект множественных сравнений (multiple comparisons/testing).
Рассмотрим это на примере. Предположим, что у нас есть 3 группы (A, B и С), в которых мы хотим сравнить среднее значение переменной интереса. Как и ранее, мы будем использовать t-тест Стьюдента. Если мы получили достаточно большое значение t-статистики такое, что p-value < 0.01, то мы отклоняем нулевую гипотезу и заключаем, что группы статистически различаются по переменной интереса. Отсечка p-value < 0.01 значит, что вероятность ошибочного вывода о различии между групповыми средними не превышает 0.01. Это будет работать именно так, когда у нас всего две группы, но в случае множественного тестирования вероятность будет больше 1%.
Выполняя тест Стьюдента, исследователь проверяет нулевую гипотезу об отсутствии разницы между двумя группами. Сравнивая группы A и В, он может ошибиться с вероятностью 1%, В и С – 1%, А и С – тоже 1%. Соответственно, вероятность ошибиться хотя бы в одном из этих трех сравнений составит:
\(P = 1 - \left(1-\alpha \right)^n = 1 - 0.99^3 \approx 0.03 > 0.01\)
Если бы групп было бы 5:
\(P = 1 - \left(1-\alpha \right)^n = 1 - 0.99^{10} \approx 0.1 > 0.01\)
К счастью, существует несколько методов, позволяющих преодолеть эту сложность.
Предположим, что мы проверяем \(n\) гипотез. Для каждой гипотезы мы будем проводить тест Стьюдента. Результаты наших тестов можно обобщить следующим образом:
| Число принятых гипотез | Число отвергнутых гипотез | Всего | |
|---|---|---|---|
| Число верных гипотез | \(U\) – Число безошибочно принятых гипотез (true negatives) | \(V\) – Число ошибочно отвергнутых гипотез (false positives) – ошибка первого рода | \(U+V=n^*\) – Число верных нулевых гипотез (true null hypotheses) |
| Число неверных гипотез | \(T\) – Число ошибочно принятых гипотез (false negatives) – ошибка второго рода) | \(S\) – Число безошибочно отвергнутых гипотез | \(T+S=n-n^*\) – Число истинных альтернативных гипотез (true alternative hypotheses) |
| Всего | \(U+T=W\) – Общее число принятых гипотез | \(V+S=R\) – Общее число отвергнутых гипотез |
В теории всего существует \(n^*\) верных нулевых гипотез. В результате наших тестов мы ошибочно отвергаем \(V\) гипотез и верно принимаем остальные \(U\) гипотез. Также существует \(n−n^*\) альтернативных гипотез, из которых \(S\) гипотез безошибочно отвергаются, а \(T\) гипотез – ошибочно принимаются. Важно, что общие количества отвергнутых и принятых гипотез (\(R\) и \(W\)), а следовательно, и суммарное число гипотез \(n\) нам известны, тогда как остальные величины (\(n^*\), \(U\), \(V\), \(T\) и \(S\)) мы не наблюдаем.
При одновременной проверке семейства статистических гипотез мы хотим, чтобы количество наших ошибок (\(V\) и \(T\)) было минимальным. Традиционно исследователи пытаются минимизировать величину ошибочно отвергнутых гипотез \(V\). Это вполне логично, поскольку ложно отвергнутая нулевая гипотеза грозит нам ложноположительным найденным эффектом, которого реально может не быть.
Если \(V \geq 1\), мы совершаем как минимум одну ошибку первого рода. Вероятность допущения такой ошибки при множественной проверке гипотез называют групповой вероятностью ошибки (familywise error rate, FWER или experiment-wise error rate). По определению, \(FWER = P(V \geq 1)\). Соответственно, когда мы говорим, что хотим контролировать групповую вероятность ошибки на определенном уровне значимости \(\alpha\), мы подразумеваем, что должно выполняться неравенство \(FWER \leq \alpha\).
Ниже мы обсудим методы, которые позволяют это делать.
Вернемся к нашему примеру, когда мы сравнили 3 группы A, B и C с помощью t-теста. Предположим, что мы получили следующие Р-значения: 0.001, 0.01 и 0.04.
Как было сказано вышем, мы хотим, чтобы групповая вероятность ошибки была не больше уровня значимости \(FWER \leq \alpha\). Согласно методу Бонферрони, мы должны сравнить каждое из полученных p-значений не с \(\alpha\), а с \(\frac{\alpha}{n}\), где \(n\) – число проверяемых гипотез. Деление исходного уровня значимости \(\alpha\) на \(n\) – это и есть поправка Бонферрони. В рассматриваемом примере каждое из полученных p-значений необходимо было бы сравнить с \(\frac{\alpha}{n}\), например, с \(\frac{0.01}{3}\approx 0.0033\).
Вместо деления уровня значимости на число гипотез, мы могли бы умножить каждое p-значение на это число и получить точно такие же выводы (эта эквивалентная процедура реалирована в R):
Иногда при домножении p-значений результат может получиться больше единицы. Из теории вероятностей мы знаем, что вероятность не может быть больше одного, поэтому в таких случаях p-значение принимают равным за единицу.
Различные виды коррекций p-значений представлены в функции p.adjust(), выбрать тип коррекции можно с помощью аргумента method. В этой функции используется домножение исходных p-значений на количество тестируемых гипотез, а не корректировка уровня значимости.
Проверим наши рассчеты:
p.adjust(c(0.001, 0.01, 0.04), method = "bonferroni")
## [1] 0.003 0.030 0.120
Можно на выходе сразу получить выводы относительно гипотез при \(\alpha =5%\):
alpha <- 0.05
p.adjust(c(0.001, 0.01, 0.04), method = "bonferroni") < alpha # отклоняем H_0 (есть эффект)?
## [1] TRUE TRUE FALSE
Важно помнить об уязвимости коррекции Бонферрони – с ростом числа гипотез мощность метода уменьшается. Чем больше гипотез мы хотим проверить, тем сложнее нам будет их отвергать (даже если они реально должны быть отвергнуты). Например, для 5 групп (10 гипотез), применение поправки Бонферрони привело бы к снижению исходного уровня значимости до 0.01/10 = 0.001. Соответственно, для отклонения гипотез, соответствующие p-значения должны быть меньше 0.001, а это довольно жесткая отсечка. Из этого делаем вывод, что при большом числе гипотез коррекцию Бонферрони лучше не использовать.
Метод Хольма позволяет побороть недостатки метода Бонферрони. Он устроен следующим образом:
Метод Хольма называют нисходящей (step-down) процедурой. Он начинается с наименьшего p-значения в упорядоченном ряду и последовательно “спускается” вниз к более высоким значениям. На каждом шаге соответствующее значение \(p-value_i\) сравнивается со скорректированным уровнем значимости \(\displaystyle{\alpha_{adjusted}=\frac{\alpha}{n+i-1}}\). Аналогично коррекции Бонферрони можно вместо корректировки уровня значимости корректировать p-значения \(\displaystyle{p-value_{i,adjusted}=p-value_{i}\cdot(n-i+1)}\) (эта эквивалентная процедура реалирована в R). Возвращаясь к нашему примеру:
А теперь проверим себя с помощью R:
p.adjust(c(0.001, 0.01, 0.04), method = "holm")
## [1] 0.003 0.020 0.040
И результаты проверки гипотез при \(\alpha =5%\):
alpha <- 0.05
p.adjust(c(0.001, 0.01, 0.04), method = "holm") < alpha # отклоняем H_0 (есть эффект)?
## [1] TRUE TRUE TRUE
Про другие процедуры family-wise error rate можно почитать тут.
Рассмотренные выше FWER методы обеспечивают контроль над групповой вероятностью ошибки первого рода. Как мы выяснили, эти методы чересчур жестко работают, когда нужно одновременно проверить слишком много гипотез (падает статистическая мощность).Под “недостаточной мощностью” понимается сохранение многих нулевых гипотез, которые потенциально могут представлять исследовательский интерес и которые, соответственно, следовало бы отклонить. Недостаточная мощность традиционных процедур множественной проверки гипотез привела к разработке новых методов, например, метода Бенджамини-Хохберга.
Для преодоления недостаточной мощности FWER методов был предложен новый подход к проблеме множественных проверок статистических гипотез. Суть подхода заключается в том, что вместо контроля над групповой вероятностью ошибки первого рода выполняется контроль над ожидаемой долей ложных отклонений (false discovery rate, FDR) среди всех отклоненных гипотез.
В терминах таблицы выше эта ожидаемая доля может быть записана следующим образом: \(\displaystyle{FDR=\left(\frac{V}{R}\right)}\).
В отличие от уровня значимости \(\alpha\), каких-либо общепринятых значений FDR не существует. Многие исследователи по аналогии контролируют FDR на уровне 5%. Интерпретация порогового значения FDR очень проста: например, если в ходе анализа данных отклонено 1000 гипотез, то при q=0.10 ожидаемая доля ложно отклоненных гипотез не превысит 100.
В статье (Benjamini, Hochberg, 1995) описание процедуры контроля над FDR выглядит так:
Эквивалентная процедура, реалированая в R отличается тем, что вместо нахождения максимального индекса \(k\), исходные p-значения корректируются следующим образом: \(q_i=\frac{p_in}{i}\).
В качестве примера рассмотрим следующий ряд из 15 упорядоченных по возрастанию p-значений (из оригинальной статьи Benjamini and Hochberg 1995):
p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BH")
## [1] 0.00150000 0.00300000 0.00950000 0.03562500 0.06030000 0.06385714
## [7] 0.06385714 0.06450000 0.07650000 0.48600000 0.58118182 0.71487500
## [13] 0.75323077 0.81321429 1.00000000
И результаты проверки гипотез при \(\alpha =5%\):
alpha <- 0.05
p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BH") < alpha # отклоняем H_0 (есть эффект)?
## [1] TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE
Интерпретация этих Р-значений с поправкой (в большинстве литературных источников их называют q-значениями) такова:
Коррекция Р-значений по методу Беньямини-Хохберга работает особенно хорошо в ситуациях, когда необходимо принять общее решение по какому-либо вопросу при наличии информации (=проверенных гипотез) по многим параметрам.
Следует помнить, что описанный здесь метод контроля над ожидаемой долей ложных отклонений предполагает, что все тесты, при помощи которых получают p-значения, независимы. На практике в большинстве случаев это условие выполняться не будет.
Для преодоления ограничения независимости тестов при проверке гипотез в работе (Benjamini and Yekutieli 2001) был предложен усовершенствованный метод, учитывающий наличие корреляции между проверяемыми гипотезами.
Процедура Бенджамини-Йекутили очень похожа на процедуру Бенджамини-Хохберга. Основное отличие заключается во введении поправочной константы \(\displaystyle{c_n=\sum \limits_{i=1}^{n}\frac{1}{i}}\), далее аналогично:
В R реализуется эквивалентная процедура:
Эквивалентная процедура, реалированая в R отличается тем, что вместо нахождения максимального индекса \(k\), исходные p-значения корректируются следующим образом: \(\displaystyle{q_i=\frac{p_i\cdot n\cdot c_n}{i}}\).
p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BY")
## [1] 0.004977343 0.009954687 0.031523175 0.118211908 0.200089208 0.211892623
## [7] 0.211892623 0.214025770 0.253844518 1.000000000 1.000000000 1.000000000
## [13] 1.000000000 1.000000000 1.000000000
И результаты проверки гипотез при \(\alpha =5%\):
alpha <- 0.05
p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BY") < alpha # отклоняем H_0 (есть эффект)?
## [1] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE
Сравним как работают разные методы:
alpha <- 0.05
n <- 50
set.seed(123)
x <- rnorm(n, mean = c(rep(0, n/2), rep(3, n/2)))
pval <- round(2*pnorm(sort(-abs(x))), 3)
default_bool <- pval < alpha
bonferroni_pval <- p.adjust(pval, method = "bonferroni")
bonferroni_bool <- p.adjust(pval, method = "bonferroni") < alpha # отклоняем H_0 (есть эффект)?
holm_pval <- p.adjust(pval, method = "holm")
holm_bool <- p.adjust(pval, method = "holm") < alpha # отклоняем H_0 (есть эффект)?
bh_pval <- p.adjust(pval, method = "BH")
bh_bool <- p.adjust(pval, method = "BH") < alpha # отклоняем H_0 (есть эффект)?
by_pval <- p.adjust(pval, method = "BY")
by_bool <- p.adjust(pval, method = "BY") < alpha # отклоняем H_0 (есть эффект)?
methods <- cbind(default_bool, bonferroni_bool, holm_bool, bh_bool, by_bool)
colnames(methods) <- c('Без коррекции', 'Бонферрони', 'Хольм', 'Бенджамини-Хохберг', 'Бенджамини-Йекутили')
rownames(methods) <- c(1:n)
methods
## Без коррекции Бонферрони Хольм Бенджамини-Хохберг Бенджамини-Йекутили
## 1 TRUE TRUE TRUE TRUE TRUE
## 2 TRUE TRUE TRUE TRUE TRUE
## 3 TRUE TRUE TRUE TRUE TRUE
## 4 TRUE TRUE TRUE TRUE TRUE
## 5 TRUE TRUE TRUE TRUE TRUE
## 6 TRUE TRUE TRUE TRUE TRUE
## 7 TRUE TRUE TRUE TRUE TRUE
## 8 TRUE TRUE TRUE TRUE TRUE
## 9 TRUE TRUE TRUE TRUE TRUE
## 10 TRUE TRUE TRUE TRUE TRUE
## 11 TRUE FALSE TRUE TRUE TRUE
## 12 TRUE FALSE FALSE TRUE TRUE
## 13 TRUE FALSE FALSE TRUE FALSE
## 14 TRUE FALSE FALSE TRUE FALSE
## 15 TRUE FALSE FALSE TRUE FALSE
## 16 TRUE FALSE FALSE TRUE FALSE
## 17 TRUE FALSE FALSE TRUE FALSE
## 18 TRUE FALSE FALSE TRUE FALSE
## 19 TRUE FALSE FALSE TRUE FALSE
## 20 TRUE FALSE FALSE TRUE FALSE
## 21 TRUE FALSE FALSE FALSE FALSE
## 22 TRUE FALSE FALSE FALSE FALSE
## 23 FALSE FALSE FALSE FALSE FALSE
## 24 FALSE FALSE FALSE FALSE FALSE
## 25 FALSE FALSE FALSE FALSE FALSE
## 26 FALSE FALSE FALSE FALSE FALSE
## 27 FALSE FALSE FALSE FALSE FALSE
## 28 FALSE FALSE FALSE FALSE FALSE
## 29 FALSE FALSE FALSE FALSE FALSE
## 30 FALSE FALSE FALSE FALSE FALSE
## 31 FALSE FALSE FALSE FALSE FALSE
## 32 FALSE FALSE FALSE FALSE FALSE
## 33 FALSE FALSE FALSE FALSE FALSE
## 34 FALSE FALSE FALSE FALSE FALSE
## 35 FALSE FALSE FALSE FALSE FALSE
## 36 FALSE FALSE FALSE FALSE FALSE
## 37 FALSE FALSE FALSE FALSE FALSE
## 38 FALSE FALSE FALSE FALSE FALSE
## 39 FALSE FALSE FALSE FALSE FALSE
## 40 FALSE FALSE FALSE FALSE FALSE
## 41 FALSE FALSE FALSE FALSE FALSE
## 42 FALSE FALSE FALSE FALSE FALSE
## 43 FALSE FALSE FALSE FALSE FALSE
## 44 FALSE FALSE FALSE FALSE FALSE
## 45 FALSE FALSE FALSE FALSE FALSE
## 46 FALSE FALSE FALSE FALSE FALSE
## 47 FALSE FALSE FALSE FALSE FALSE
## 48 FALSE FALSE FALSE FALSE FALSE
## 49 FALSE FALSE FALSE FALSE FALSE
## 50 FALSE FALSE FALSE FALSE FALSE
plot(pval, bonferroni_pval, col = "orange", type="p", pch=1)
lines(pval, holm_pval, col="green", type="p", pch=1)
lines(pval, bh_pval, col="blue", type="p", pch=1)
lines(pval, by_pval, col="violet", type="p", pch=1)
abline(h=alpha, col="red")
abline(v=alpha, col="red")
legend(x=0.6, y=0.5,
legend=c('Бонферрони', 'Хольм', 'Бенджамини-Хохберг', 'Бенджамини-Йекутили', 'Уровень значимости'),
col=c("orange", "green", "blue", "violet", "red"),
bty = "n", pch=1)
Ранее мы говорили, что при достаточном объеме выборки и хорошей рандомизации в среднем тритмент группа и контрольная группа окажутся похожимы друг на друга за исключением оказываемого на них воздействия, в этом случае в качестве хорошей оценки эффекта можно использовать простую разницу в средних исходах.
По разным причинам может получиться так, что тритмент группа и контрольная группа не похожи друг на друга.
Существует 2 способа как бороться с несбалансированностью групп:
Если у вас возникли какие-то вопросы, или вы нашли неточности в тексте, напишите мне об этом на почту annastavnychuk@gmail.com